function [adaptationCurve] = stimHistory(subjID, t_thresh)

% PS: If t_thresh is greater than 10, it is intepreted as the desired
% number of voxels that the ROI should have, and the t-value is calculated

adaptationCurve = cell(4,7);

adaptationCurve{1,1} = subjID;
adaptationCurve{1,2} = '1-back';
adaptationCurve{1,3} = '2-back';
adaptationCurve{1,4} = '3-back';
adaptationCurve{1,5} = '4-back';
adaptationCurve{1,6} = '5-back';
adaptationCurve{1,7} = '6-back';

adaptationCurve{2,1} = 'FOBS mask';
adaptationCurve{3,1} = 'Adapt mask';
adaptationCurve{4,1} = 'Intersection';


% Define some basic directories and files
subjectDir = ['/Volumes/cluster/jet/mattar/AdaptID/Analyses/' subjID '/'];
FOBSdir = [subjectDir 'FOBSLoc.feat/'];
FOBSfileZIPED = [FOBSdir 'stats/tstat1.nii.gz'];

% Gunzip FOBS data
FOBSfile_filename = gunzip(FOBSfileZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
FOBSfile_filename = FOBSfile_filename{1};

% Load FOBS data .nii file in MATLAB
FOBSroiStruct = load_untouch_nii(FOBSfile_filename);
FOBSroi = FOBSroiStruct.img;
if t_thresh > 10
    [~,sortIndex] = sort(FOBSroi(:),'descend');
    FOBSroi_mask = zeros(size(FOBSroi));
    FOBSroi_mask(sortIndex(1:t_thresh)) = 1;
else
    FOBSroi_mask = FOBSroi > t_thresh;
end

maskSize_FOBSroi = sum(FOBSroi_mask(:));
cubeDims = size(FOBSroi_mask);

% Save FFA mask in subject's masks directory
FFAmaskStruct = FOBSroiStruct;
FFAmaskStruct.img = FOBSroi_mask;
if exist([subjectDir 'masks'],'dir') == 0
    mkdir([subjectDir 'masks']);
end
save_untouch_nii(FFAmaskStruct, [subjectDir 'masks/FFAmask.nii']);

% Delete FOBS .nii file
eval(['delete ' FOBSfile_filename]);



% Retrieve the analysis directories
FaceDirectories = struct2cell(dir([subjectDir 'Face*+.feat']));
FaceDirectories = FaceDirectories(1,:)';
numRuns = length(FaceDirectories);

% Allocate variables to save all relevant maps
cope_feat = cell(numRuns,6);
varcope_feat = cell(numRuns,6);
tstat_feat = cell(numRuns,6);

cope_gfeat = cell(1,6);
varcope_gfeat = cell(1,6);
tstat_gfeat = cell(1,6);

% Retrieve cope, varcope and tstat maps
for runIndx = 1:numRuns
    for stepsBack = 1:6
        % Generate filenames
        copeZIPED = [subjectDir FaceDirectories{runIndx} '/stats/cope' num2str(stepsBack+1) '.nii.gz'];
        varcopeZIPED = [subjectDir FaceDirectories{runIndx} '/stats/varcope' num2str(stepsBack+1) '.nii.gz'];
        tstatZIPED = [subjectDir FaceDirectories{runIndx} '/stats/tstat' num2str(stepsBack+1) '.nii.gz'];
        
        % Gunzip cope, varcope and tstat
        cope_filename = gunzip(copeZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
        cope_filename = cope_filename{1};
        varcope_filename = gunzip(varcopeZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
        varcope_filename = varcope_filename{1};
        tstat_filename = gunzip(tstatZIPED, '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis');
        tstat_filename = tstat_filename{1};
        
        % Load data into matlab
        this_cope = load_untouch_nii(cope_filename);
        this_cope = this_cope.img;
        this_varcope = load_untouch_nii(varcope_filename);
        this_varcope = this_varcope.img;
        this_tstat = load_untouch_nii(tstat_filename);
        this_tstat = this_tstat.img;
        
        % Save maps
        cope_feat{runIndx,stepsBack} = this_cope;
        varcope_feat{runIndx,stepsBack} = this_varcope;
        tstat_feat{runIndx,stepsBack} = this_tstat;
        
        % Delete gunziped files
        eval(['delete ' cope_filename]);
        eval(['delete ' varcope_filename]);
        eval(['delete ' tstat_filename]);
    end
end


% Generate the t-stat map for the group
for stepsBack = 1:6
    sumofoneovervarcope = zeros(cubeDims);
    sumofcopeovervarcope = zeros(cubeDims);
    for thisRun = 1:numRuns
        sumofoneovervarcope = sumofoneovervarcope + 1./varcope_feat{thisRun,stepsBack};
        sumofcopeovervarcope = sumofcopeovervarcope + cope_feat{thisRun,stepsBack}./varcope_feat{thisRun,stepsBack};
    end
    varcope_gfeat{1,stepsBack} = 1./sumofoneovervarcope;
    cope_gfeat{1,stepsBack} = (1./sumofoneovervarcope) .* sumofcopeovervarcope;
    cope_gfeat{1,stepsBack}(isnan(cope_gfeat{1,stepsBack})) = 0;
    tstat_gfeat{1,stepsBack} = cope_gfeat{1,stepsBack} ./ sqrt(varcope_gfeat{1,stepsBack});
    tstat_gfeat{1,stepsBack}(isnan(tstat_gfeat{1,stepsBack})) = 0;
end

   


    


% Generate the mask for calculating average effect
if t_thresh > 10
    [~,sortIndex] = sort(tstat_gfeat{1,1}(:),'descend');
    FaceAdapt_mask = zeros(cubeDims);
    FaceAdapt_mask(sortIndex(1:t_thresh)) = 1;
else
    FaceAdapt_mask = tstat_gfeat{1,1} > t_thresh;
end
maskSize_FaceAdapt = sum(FaceAdapt_mask(:));
    
    
% Now, calculate the effects on the intersection of the matrices
% Generate the mask for calculating average effect
minFFAMatrix = min(FOBSroi,tstat_gfeat{1,1});
if t_thresh > 10
    [~,sortIndex] = sort(minFFAMatrix(:),'descend');
    intersectionFFA_mask = zeros(cubeDims);
    intersectionFFA_mask(sortIndex(1:t_thresh)) = 1;
    resulting_t_threshold = minFFAMatrix(sortIndex(t_thresh));
else
    intersectionFFA_mask = minFFAMatrix > t_thresh;
end
maskSize_intersectionFFA = sum(intersectionFFA_mask(:));




% Now, calculate the average betas for the stimulus history covariates
for stepsBack = 1:6
    adaptationCurve{2,stepsBack+1} = sum(cope_gfeat{1,stepsBack}(:) .* FOBSroi_mask(:)) / maskSize_FOBSroi;
    adaptationCurve{3,stepsBack+1} = sum(cope_gfeat{1,stepsBack}(:) .* FaceAdapt_mask(:)) / maskSize_FaceAdapt;
    adaptationCurve{4,stepsBack+1} = sum(cope_gfeat{1,stepsBack}(:) .* intersectionFFA_mask(:)) / maskSize_intersectionFFA;
end

save(['workspace_' subjID '.mat']);
